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ABSTRACT 

It has been previously shown that the temporal integration of hyperbolic partial differ- 
ential equations may , because of boundary conditions, lead to deterioration of accuracy of 
the solution. A procedure for removal of this error in the linear case has been established 
previously. 

In the present paper we consider hyperbolic p.d.e’s (linear and non-linear) whose bound- 
ary treatment is done via the SAT-procedure. A methodology is present for recovery of the 
full order of accuracy, and has been applied to the case of a 4th order explicit finite difference 
scheme. 
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1. Introduction 

Currently there is a growing interest in long time integration for solving problems in areas 
such a aero-acoustics, electro-magnetics, material-science and others. This necessitates, ([1]), 
working with higher order (4th order accuracy and above) schemes. Often the methodology 
of choice is to semi-discretize the equations by applying a high order (4th and above) spatial 
difference operator and then advance temporally using single level multi-stage Runge-Kutta 
integrators. This raises the question of how to supply boundary values at the intermediate 
stages of the R-K integrations. In the case of hyperbolic p.d.e’s that means the imposition 
of the time dependent conditions at the inflow boundary. 

The conventional (and intuitively natural) way of imposing inflow boundary conditions 
at the intermediate stages is to use the “appropriate” value of the boundary data, g(t), at 
each stage. Thus, for example, at a stage corresponding to t + (A</2), one would impose 
g(t + At/2). 

In a previous paper ([2]) , it was shown that the procedure described above, when ap- 
plied to hyperbolic p.d.e’s with time dependent b.c.’s, reduces the accuracy near the inflow 
boundary to first order and thus the overall accuracy cannot exceed 0(A.t 2 ). This conclusion 
is independent of the order of accuracy of the spatial difference operator. 

One way of avoiding the dilemma of what boundary values one should supply at the 
intermediate stages, is to advance the R-K integration without imposing any intermediate 
values, but rather obtain the intermediate boundary values from the numerical solution 
operator. However, this approach has the disadvantage of reducing substantially the stability 
limit (e.g., the allowable time step is reduced by a factor 2 in the case of 4th order classical 
R-K with a 4th order spatial derivative operator), hence rendering it less than attractive. 

In ([2]) a general methodology was presented in the case of linear p.d.e.’s for the correct 
imposition of the intermediate stage boundary values so that the scheme recovers its full 
formal accuracy. This was expounded in detail for the case of the classic 4th order R-K 
with a hyperbolic 4th order spatial difference operator. It was also shown there that in the 
non-linear case (e.g. hyperbolic conservation laws) this methodology was applicable to R-K 
integration up to 3rd order . For R-K . methods of 4th order and above we were not able to 
extend the theoretical approach described in ( [2]). 

In this paper we address anew the issue of how to deal with the non-linear case. We 
present a methodology for retaining the full accuracy even in the non-linear case. The 
application of this methodology involved numerical determination of free parameters in con- 
tradistinctsr numerical determination of free parameters in contradistinction to the linear 
procedure described in ([2]). We find for example, that in the 4th order classical R-K integra- 
tor with 4th order explicit spatial derivative operator, the full accuracy is retained without 
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any reduction in t lie allowable time step. 

The new procedure is demonstrated for hyperbolic problems where the boundary condi- 
tions are satisfied by using the SAT approach ([3]. The motivation for doing so is that the 
SAT procedure is the only one that prevents temporal growth not present in the true solu- 
tion of a system of p.d.e.’s. Section 2 describes how to apply correctly the intermediate SAT 
boundary conditions in the case of a linear problem. In Section 3 we cover the non-linear 
case. 


2. The Linear Case 

In this section we analyze the effect of imposing the inflow boundary conditions in the 
conventional way when the discretization algorithm employs the SAT approach (see ([3]) ). 
The SAT is a penalty type method that was constructed so as to ensure that the numerical 
solution will not include temporal growth which is not of a physical origin. This is achieved 
by mimicking the energy estimate of the p.d.e. 

Recall that we are considering the following hyperbolic problem (see ([3])): 


du du 

a7 + & = °- 

(2.1) 

II 

o 

<-> 

(2.2). 


The SAT formulation for the semi-discrete version of (2.1) - (2.2), based on a uniform grid, 
is 

^ = (iflj7(0).-ir?Mf)-sW). >' = 0,1 N; t> 0 (2.3) 

where V= [t’oi in, . . . , ujv] T is the semi-discrete approximation that converges to n(xi,t) at 
the spatial grid points (for stable discretizations); and is the differential matrix rep- 
resentation of the derivative operator (— £;). The vector <r depends on the differentiation 
matrix j^D, and on the energy norm used in bounding the error. It is determined as described 
in ([3]); see the discussion after equation (6) therein. 

The demonstration of accuracy deterioration will be shown for the four stage “classical” 
RK algorithm, which is one of the most commonly used RK time advancing schemes. For 
the analysis to make sense we assume that the spatial discretization is at least fourth order 

I 

accurate. 

The above mentioned four stage RK integrator is implemented as follows: 

I /(1) = y(") + ± DV (n) _ A r ^ [u< n) - g(t ) ] (2.4) 
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(2.5) 


V m = K'"» + ±DV' "1 - JT (7 [4'> - g (t + y)] 

(/(3| = ,/W + - At 5 t>J 2) - g (t + y)l (2.G) 

J/(n+l) _ y( n) + A /)[]/(*») _|_ 2I /(1) + 2V (2) + V (3) ] 

6 

-^{h’o n) -<7(0] + 2 v { 0 l) -g(i + y)] + 2 + 

+bo 3) — 50 + AO]} • 

To check for accuracy we substitute for the V ^ the exact values rt(x,-, t), and in particular 
a= g[t). Note from eq. (2.4), that on the boundary (using the differential equations 
( 2 . 1 ), ( 2 . 2 )) 

•4" = >4"’ - y (£“) o = + y»'«- (2 ' 8) 

Thus in eq. (2.7), we have for the term — g (t + ^), 

4° - s (‘ + y) = [*M + ys'w] - s (‘ + y) = °( A ‘ 2 )- < 2 - 9 ) 

Thus y( n+1 ) — f/( n ) is at best 0(Af 2 ), and not 0(Af 5 ) as expected from the R-K scheme used. 

In this linear case, the remedy proposed in the previous paper, ([2]), works here as well. 
In particular, eqs. (2.4) -(2.7) take the following form: 


V<" = I 7 '” 1 + - j(<) ! 

(2.10) 

W = V<"> + ~DV^ - It a [4° - g(t) - ys'(f)] 

(2.11) 

K (3) = V 7 '”’ + A DV m - Ar a .-J 2 ' - ,(1) - y g'(t) - yV(«)l 

(2.12) 

yn+l _ yn + Jj/(«) + 2 |/0) + 2V (2) + V (3) ] 


-^{[r4 n) - 5(01 + 2 [vo ) ~5(0 “ y5'0)] 


+2 |t>o — 5(0 “ 2 5'(0“ 4 5 (0 

(2.13) 


o) / \ a / At 2 .. A 
+ t’o - 5(0 ~ 9 - —9 ~ —9 

It is readily verified that V r f n+1 ^ — = 0(A< 5 ), as required. 
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3. The Non-Linear Case 

For the sake of simplicity we consider first the scalar conservation law p.d.e. 

757 + ^=°' <>o (3.1) 

at ax 

«(0 ,t)=g(t). (3.2) 

In general, for any spatial discretization (whether explicit or implicit) the semi-discrete 
form of (3.1) - (3.2) is: 

= in/(V) - ir ? [,* - ,(()]. (3.3) 

Using the notation of reference ([3]), 

D = -P~'Q (3.4) 

where \D is the differentiation matrix representing the differential operator, (-Jj); com- 
posed of the explicit part Q and the inverse of the implicit part P. For a fully explicit spatial 
differentiation, P = / + J3, where B differs from zero only at the two diagonal corners (see 
examples of P and Q in ([4])).. 

The vector a, again using the notation of [ ], is given by 

a= ha{u 0 )g 0 oP- l H- l S (3.5) 

where a(u 0 ) = (f{) Q , and goo is twice the value of the left upper corner element of HQ. 
For the definition of the matrix H see Assumption I in [1]. The parameter r is determined 
from stability consideration to be r > 1, see ([3]). 

Next we demonstrate that writing the classical 4th order RK for eq. (3.3), using the 
linear “fix” as in equations (2.10) - (2.13), does not yield the required 4th order accuracy: 

1/(11 = K<»> + jD/(!/ |n| ) - |r ? (4"> - «(()] (3.5) 

= v<"> + ^/(F"’) - a [4" - ,(() - y S '(i)] (3.6) 

A 4 A. 4 2 

V<»> = F<”> + A D/(F«) - Ar 5 v™ - g(t) - “<,'(() - ^-g"(t) (3.7) 
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0n+l) = y(n) + 2-D[f(V n ) + 2f(V {1) ) + 2f(V {2) ) + f(VW)} 

~^ T(T {(4 n) - 5(0] + 2 v 0 } - 5(0 - ^•5 , (0] 


+2 uj 2) - g{t) - y5'(0 - ^p5"(0 


(3.8) 


+ 


(3) /i\ Ai ' ^ n" ^ 

I’o - 5(0 - A<5 - —5 - 


]}■ 


2 4 

Again, when checking accuracy, we take = u(.r,-, <), and in particular v ^ = g(t ); and 
also \Df(V n ) = -£/(u) + 0 (Af 4 ). With these preliminaries we get immediately from 
equation (3.5) 


« 4 " = vS - y A/W. = 9(<) + y|(“)o = 9(1) + Y9'(D- 


(3.9) 


Note that this is the same as in the linear case, see eq. (2.8). Thus with = jf(0 an d 
with equation ( 3 . 9 ) we can supply, for the purpose of accuracy checking, the correct values 
of t> (n) and v^\ When we look at eq. (3.6), using the above results, the governing non-linear 
p.d.e., and simple Taylor’s expansion, we have: 


y(2) = ]A n) + -Df 
2 


V {n) + ^Df(V n ) 


= v {n) ~— — f 

2 dx J 


„-f| /( „) + 0 ( A i , ) 


= v n - — — / 

2 dx J 


At dll n , . ,r. 

u + T¥ + 0(A() 




= u 


At 

■ 

At 1 

U ~Y 

-lit - 

- Y u tt \ 


+ 0 (Af 3 ). 
So finally we have on the boundary 

,i 2 ) = n(l) 4- —n' 


At 2 


»o = 5(0 + -7T9 \ 0 + ~r5"(0 + °( A * ) 


(3.10) 


(3.11) 


o 2 J ’ 4 

It follows from (3.11) that the “penalty” term in (3.8) introduces an error of (At 3 ). Since 
the coefficient of (At 3 ), [/ uu (u<) 2 ]o> cannot generally be expressed as a function of #(0 and 
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its derivatives, it is very difficult to remedy the situation. Thus the “linear procedure” fails 
at the third RK stage. 

We now propose a methodology to deal with the RK integration of non linear hyperbolic 
conservation laws. We first present this procedure in the case of the classic 4th order R-K 
scheme. Our starting point is the observation that the “linear-procedure” yields the required 
accuracy for i>q and Vq\ The idea is to use at each stage a linear combination of the “linear” 
SAT, or penalty, terms used (3.5) and (3.6). The 4th order classic RK stages will thus be: 

V (1) = yM + Ad/(V) - a [ti n) - flr(O) (3.12) 

V (2) = y (n) + ^Df{ 1 /(1) ) - ^ft a [4° - g(t) - Y 9 '^ ~ p * f u ° n) “ $(*)! ( 3 - 13 ) 

y(3) = V (n) + A Df(V { V) - A 6 a [i4 n) - g(t)) - Ae a [uj 1 * - g(t) - y </'(*)] (3.H) 

^(n+i) = vM + ^D[f{V n ) + 2f{VW) + 2f(VW) + f(VW)} 

(3.15) 

-pt V [vo ] - 5(0] - p ° } - g(t) - y5(0 

where the free parameters a, ft, 7, 6, e , ft, and u will be chosen so as to maximize the allowable 
time step. It is clear from the previous discussion that the system (3.12) - (3.15) maintains 
the 4th order accuracy. There remains the question of whether the CFL stability condition 
deteriorates, in comparison to the conventional application of R-K, equation (2.4) - (2.7). It 
is also clear that the optimal choice of the free parameters varies with the spatial 

discretizations (i.e. the differentiation matrix D) and boundary closures. One checks that 
the absolute values of the eigenvalues of the amplification matrix resulting from (3.12) - 
(3.15) should not exceed unity. We carried out this procedure (using Matematica Software) 
in the case of an explicit 4th order algorithm with 3rd order boundary closures. In this case 
// = /, and the matrices P and Q can be found in section (9.1) of reference [3]. With the 
following values of the free parameters, a — ft = —e = 1, 6 — 2, ft = 0, v = 3, 7 = —0.37, 
the CFL condition becomes A < 2.1. This is the same restriction on the time step, At, as 
one has in the linear problem using the “conventional”, i.e. less accurate, boundary values 
with or without the SAT-term. 

4. Conclusions 

In summary, we have in (3.12)-(3.15) a 4th order RK scheme, applied to a non-linear 
p.d.e., which maintains the overall 4th order accuracy without any decrease in the allowable 
time step. The extension to a system of hyperbolic p.d.e.’s is quite straight forward using 
the SAT-system approach delineated in ([3]). 
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